In [1]:
% pylab inline
from scipy.optimize import curve_fit
Fermi-Dirac distribution, g=2¶
In [2]:
# The main function for determining the microstates in a Fermi system with possible degeneracy
# but it represents the states listing the energy levels of the particles
# this function is used by 'occupationsFg'
def disg(n,N,g=1,m=-1,nm=0,d=0):
'''determines the possibilities,
how n excitations can be distributed among N particles
in states with energies k+1/2, k=0,1,2,...,
and maximum g number of particles in each state,
and gives out a list of the energies of the particles for each possibility.
For the recursive use of the function the lowest energy level index m with
already fixed particle[s] and its occupation number nm can be given.'''
# q=' '.join([' ']*d) # depth indent in test mode
s=(N-1)//g # starting level index of the particle with largest energy
ns=(N-1)%g +1 # starting occupation of that level
if m==-1: m=s+n # if m is not given choose a high enough value
# print(q,'N,n,g,m,nm,s,ns ',N,n,g,m,nm,s,ns)
if N==1:
if n<m or (n==m and nm<g): return [[n]]
else: return []
else:
f0=s # lowest possible final level
if nm<g: f1=min(m,s+n)
else: f1=min(m-1,s+n) # highest possible final level
# print(q,'f0,f1',f0,f1)
p=[]
for f in range(f0,f1+1):
# print(q,'f',f)
if f<m: p1=disg(n-(f-s),N-1,g,f,1,d+1)
else: p1=disg(n-(f-s),N-1,g,f,nm+1,d+1)
# print(q,'p1',p1)
for i in range(len(p1)):
p1[i].append(f)
# print(q,'p1[i]',p1[i])
p+=p1
return(p)
In [3]:
# The main function for determining the microstates in a Bose system
# It represents the states listing the occupation numbers of the energy levels
def occupationsFg(E,N,g=1,l=0):
'''Determines the possibilities how N particles can be positioned
in states with energies k+1/2, k=0,1,2,...,
so that the full energy is E,
and gives out a list of the occupation numbers for each possibility.
l is the length of the lists,
which is set to E-N/2+1 if not given in the function call.'''
s=(N-1)//g # index of the highest occupied state in the ground state
ns=(N-1)%g +1 # occupation of that level in the ground state
E0=g*s**2/2+ns*(s+1/2) # energy in the ground state
# print('E0 ',E0)
ne=E-E0 # the sum of excitations
n=round(ne)
if n!=ne: return [] # ne should be an integer, i.e.
# Since one particle can get at most the full excitation energy n,
# the highest possibly occupied state will be the one indexed by s+n.
# So the necessary length of the occupation number list is s+n+1.
# This shall be the default value of l
# and a given value of l lower than s+n+1 worth a warning
if l==0: l=s+n+1
if l<s+n+1: print('warning: l is too low to show all occupied states')
# the excitation energy shall be distributed among the N particles as n quantums
da=disg(n,N,g)
# 'da' is a list of possible distribution of the energy among the particles,
# each distribution is represented by a list 'di' that contains the excitations of every particle;
# for each distribution these excitations should be counted in the corresponding
# elements of the list 'oc' of occupation values
# and form a list of these 'oc' lists for all distributions
ocs=[]
for di in da:
oc=[0]*l
for i in di:
if i<l: oc[i]+=1
ocs+=[oc]
return ocs
In [4]:
N=20; E=106; g=2; l=20
nilists=occupationsFg(E,N,g,l)
nilists.reverse()
nilists
Out[4]:
In [5]:
print('N: ',N)
print('E: ',E)
s=(N-1)//g; ns=(N-1)%g +1; E0=g*s**2/2+ns*(s+1/2) # energy in the ground state
print('excitation: ',E-E0)
ei=array(range(l))+0.5
print()
#print('not taking into account spin-multiplicity')
print('number of possible configurations:',len(nilists))
#niav=sum(array(nilists),axis=0)/len(nilists)
#print('average occupation numbers:')
#print(niav)
#print('check N:',sum(niav),', E:',sum(ei*niav))
print()
print('Taking into account spin-multiplicity:')
import scipy.special
nm=0
nisum=array(nilists[0])*0
for nilist in nilists:
m=int(prod(scipy.special.binom(g,nilist)))
nm+=m
nisum+=array(nilist)*m
print('number of microstates:',nm)
niav=nisum/nm/g
print('average occupation numbers:')
print(niav)
print()
In [6]:
def Fermi(x,mu,kT):
return 1/(exp((x-mu)/kT)+1)
parameter, covariance_matrix = curve_fit(Fermi, ei, niav)
fit=Fermi(ei,*parameter)
print('fitted mu:',parameter[0],', kT:',parameter[1])
print('sum of squared deviations in fit:',sum((fit-niav)**2))
print('approximation from the fitted data for N:',sum(g*fit),', for E:',sum(ei*g*fit))
xvec=linspace(ei[0],ei[-1],100)
fitvec=Fermi(xvec,*parameter)
plot(ei,niav,"o",xvec,fitvec);
In [ ]:
Bose-Einstein distribution¶
In [7]:
# The main function for determining the microstates
# both in a Bose system and in a nondegenerate Fermi system,
# but it represents the states listing the energy levels of the particles,
# this function is used by 'occupationsB' and 'occupationsFg1'
def disall(n,c=0,m=0):
'''determines the possibilities,
how n quantums can be distributed in c columns
in a monotonically increasing way
with maximum m quantums in each column
and also prints the empty columns;
c and m are set to n if not given in function call'''
if c==0: c=n
if m==0: m=n
if c==1 or n<=1:
if n<=m:
return [[0]*(c-1)+[n]]
else:
return []
else:
p=[]
for l in range(1,min(m,n)+1):
if l<n:
p1=disall(n-l,c-1,l)
for i in range(len(p1)):
p1[i].append(l)
p+=p1
else:
p+=[[0]*(c-1)+[n]]
return(p)
In [8]:
# The main function for determining the microstates in a Bose system
# It represents the states listing the occupation numbers of the energy levels
def occupationsB(E,N,l=0):
'''Determines the possibilities how N particles can be positioned
in states with energies k+1/2, k=0,1,2,...,
so that the full energy is E,
and gives out a list of the occupation numbers for each possibility.
l is the length of the lists,
which is set to E-N/2+1 if not given in the function call.'''
# Since in the ground state E_0=N/2 the sum of excitations should be
ne=E-N/2
# ne should be an integer, i.e.
# if N is even E should be integer, if N is odd E should be half integer:
n=round(ne)
if n!=ne: return []
# Since one particle can get at most the full excitation energy n,
# the highest possibly occupied state will be the one indexed by n.
# So the necessary length of the occupation number list is n+1.
# This shall be the default value of l
# and a given value of l lower than n+1 worth a warning
if l==0: l=n+1
if l<n+1: print('warning: l is too low to show all occupied states')
# the excitation energy shall be distributed among the N particles as n quantums
da=disall(n,N)
# 'da' is a list of microstates,
# each represented by a list 'di' that contains the excitations of every particle;
# for each microstates these excitations should be counted in the corresponding
# elements of the list 'oc' of occupation values of the states
# and form a list of these 'oc' lists for all microstates
ocs=[]
for di in da:
oc=[0]*l
for i in di:
if i<l: oc[i]+=1
ocs+=[oc]
return ocs
In [9]:
N=10; E=20; l=18
nilists=occupationsB(E,N,l)
nilists
Out[9]:
In [10]:
print('N: ',N)
print('E: ',E)
print('excitation: ',E-N/2)
print('number of microstates:',len(nilists))
niav=sum(array(nilists),axis=0)/len(nilists)
print('average occupation numbers:')
print(niav)
ei=array(range(l))+0.5
print('check N:',sum(niav),', E:',sum(ei*niav))
In [11]:
def Bose(x,mu,kT):
return 1/(exp((x-mu)/kT)-1)
parameter, covariance_matrix = curve_fit(Bose, ei, niav, p0=(0,1))
fit=Bose(ei,*parameter)
print('fitted mu:',parameter[0],', kT:',parameter[1])
print('sum of squared deviations in fit:',sum((fit-niav)**2))
print('approximation from the fitted data for N:',sum(fit),', for E:',sum(ei*fit))
xvec=linspace(ei[0],ei[-1],100)
fitvec=Bose(xvec,*parameter)
plot(ei,niav,"o",xvec,fitvec);
In [ ]: